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Abstract : We study a 2D quasi-static discrete crack anti-plane model of a 
tectonic plate with long range elastic forces and quenched disorder. The plate is 
driven at its border and the load is transfered to all elements through elastic forces. 
This model can be considered as belonging to the class of self-organized models which 
may exhibit spontaneous criticality, with four additional ingredients compared to 
sandpile models, namely quenched disorder, boundary driving, long range forces and 
fast time crack rules. In this "crack" model, as in the "dislocation" version previously 
studied, we find that the occurrence of repeated earthquakes organizes the activity 
on well-defined fault-like structures. In contrast with the " dislocation" model, after a 
transient, the time evolution becomes periodic with run-aways ending each cycle. This 
stems from the " crack" stress transfer rule preventing criticality to organize in favor of 
cyclic behavior. For sufficiently large disorder and weak stress drop, these large events 
are preceded by a complex space-time history of foreshock activity, characterized by a 
Gutenberg-Richter power law distribution with universal exponent B = 1±0.05. This 
is similar to a power law distribution of small nucleating droplets before the nucleation 
of the macroscopic phase in a first-order phase transition. For large disorder and 
large stress drop, and for certain specific initial disorder configurations, the stress 
field becomes frustrated in fast time : out-of-plane deformations (thrust and normal 
faulting) and/or a genuine dynamics must be introduced to resolve this frustration. 



Resume : Nous etudions un modele quasi-statique discret anti-plan de fissures 
dans une plaque tectonique en presence de forces elastiques a longue portee et de 
desordre gele. La plaque est soumise sur ces bords a un taux de cisaillement con- 
stant faible qui est retransmis sur tons les elements de la plaque par les interactions 
elastiques. Ce modele pent etre vu comme un exemple de modeles auto-organises 
qui peuvent donner lieu a une criticalite spontanee, possedant les ingredients addi- 
tionnels suivants : desordre gele, forgage de frontiere, transfert a longues distances 
et regies d'evolution en "temps rapide" de type "fissure". Dans ce modele de type 
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"fissure", comme dans la version de type "dislocation" precedamment etudiee, les 
tremblements de terre repetcs s'organisent sur dcs structures dc failles bien definies. 
A la difference du modele de type "dislocation", on observe qu'aprcs un transitoire 
revolution temporelle devient periodique, avec des cvcnements catastrophiques qui 
terminent chaque cycle. Cela provient du transfert des contraintes a grandes dis- 
tances creees par les fissures, qui empeche la criticalite de s'organiser en faveur d'un 
comportement cyclique. Pour des desordres suffisamment grands et des baisses de 
contraintes associees aux ruptures suffisamment faibles, les grands tremblements de 
terre sont precedes d'une histoire spatio-temporelle complexe decrivant I'activite des 
precurseurs, caracterisee par une loi de puissance de type Gutenberg-Richter avec un 
exposant universel 5 = 1±0.05. Ce resultat est similaire a une distribution de petites 
"gouttes" de nucleation preccdant la formation d'une phase macroscopique dans une 
transition du premier ordre. Pour des forts desordres et de baisses de contraintes im- 
portantes lors des ruptures, et pour des realisations specifiques du desordre, le champ 
de contrainte devient "frustre" pendant un tremblement de terre en "temps rapide" 
: des deformations hors-plan (failles inverses et normales) et/ou une vrai dynamique 
doivent etre introduit pour enlever cette frustration. 



PACS : 91.30.Dk : Seismicity : space and time distribution 91.45.-c : Physics 
of plate tectonics 64.60.Ht : Dynamical critical phenomena 05.40. : Fluctuation 

phenomena, random processes and Brownian motion 
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1 Introduction 



In this paper, we continue our exploration of the hypothesis according to which 
earthquake and fault characteristics can be understood, at time scales of years and 
above, only by using a global perspective, treating on the same level the growth of 
faults by repeated earthquakes on one hand and the localization of earthquakes on 
faults on the other hand. A lot of studies have documented the self-similar structure 
of fault patterns 0, ||, |^, |, ||, |, ^ . At short time scales (typically less than or of 
the order of a century), earthquakes occur on these pre-existing set of faults and one 
can neglect the evolution of the fault network to focus on the question of the role of 
the fault structure in the observed earthquake phenomenology. On the other hand, 
faults are evolving, nucleating, growing, branching, healing, dying eventually, being 
screened by others faults ^ . This evolution occurs as a result of the accumulation of 
deformations, accounted for by earthquakes for a significant fraction (varying upon 
the location on earth). Therefore, the fault geometry, which must be introduced 
to get a correct description of earthquake occurrence, is not an arbitrary fractal, 
but results from the accumulation of earthquakes and other more ductile modes of 
deformations which are themselves determined by the geometry. Our purpose here 
is to explore further the implications of this "hen and egg" problem, within a simple 
model. 

To tackle this question, we have previously introduced a 2-D quasi-static "dislo- 
cation" model for the generation and organization of faults by repeated earthquakes 



in an heterogeneous elastic plate driven at its border |Ty, |Tl|, |I2[. The main results 
are the spontaneous generation of fractal fault structures and the existence of a well- 
defined Gutenberg-Richter earthquake energy distribution N{E)dE ~ E-^^+^^dE, 
with B = 0.3 ± 0.05 in 2D describing the earthquake population. The faults have 
been found to be globally optimal structures in the sense that they can be mapped 
onto a minimal interface problem, which in 2D corresponds to the random directed 



polymer problem [jT3|. More precisely, for small stress drops, we have shown that a 
fault minimizes the sum of the random thresholds of the elements along it. This global 
minimization problem is achieved by the spontaneous organization of the medium, 
in which after a long "learning" transient regime, the deformation becomes localized 
on the optimal fault structure. In a sense, the elastic plate, endowed with its rules of 
rupture and stress redistribution, can be viewed as an analog computer which solves 
an optimization problem. Concretely, the outcome of this optimization correspon- 
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dence is that the faults are self-affine with a roughness exponent which is known 
exactly and equal to 2/3. 

In the present work, we study the "crack" version of the model : a constant "dy- 
namical" stress is assumed to characterize the "fast time" rupture, i.e. the stress on 
all broken elements is fixed during a given earthquake event. In contrast, the previous 
"dislocation" model |T^, 0, |T2| corresponds to imposing a slip to the ruptured ele- 
ment, allowing it to be reloaded in fast time in the succession of ruptures producing 
a complete earthquake. In the present "crack" model, an element thus fails only once 
in a given earthquake and stress enhancement occurs at the crack tips, growing, as 
usual, as the square root of the length of the evolving earthquake. The motivation 
to study this variant of the initial model is twofolds: 1) there is some debate in the 
litterature on the correct model (dislocation or crack) to use for large earthquakes 
T^, |TB[ ; 2) in the context of self-organization, it is important to assess the role of 



local rules in the resulting large-scale organization fl% . 

We explore the various regimes as a function of initial disorder (on the stress 
thresholds and elastic coefficients) and dynamical stress drop amplitude. As in the 
"dislocation" model, we find that the occurrence of repeated earthquakes organizes 
the activity on well-defined fault-like structures. The main difference with the " dislo- 
cation" model stems from the tendency for the model to synchronize, i.e. to generate 
large periodic events, albeit of a complex internal spatio-temporal structure. We in- 
terpret this behavior as resulting from the physics of coupled relaxation oscillators 



with threshold [|T8|, [19], [2^, |2T], For sufficiently large disorder and weak stress 

drop, these large events are preceded by a complex space-time history of foreshock 
activity, characterized by a Gutenberg- Richter distribution with exponent B = 1±0.1 
which is universal in the sense that the exponent is found essentially the same for all 
systems explored. For large disorder and large stress drop, for certain specific initial 
disorder configuration, and after a long time, the model does not have a solution 
anymore: this corresponds to a situation where the stress field becomes "frustrated" 
and the anti-plane quasi-static crack modelling is no more defined. This shows that 
the quasi-static "crack" version of the model is not self-consistent and additional 
types of deformations must be allowed to get rid of this frustration. For instance, 
out-of-plane deformations (thrust and inverse faulting) and/or the introduction of a 
genuine dynamics can resolve this frustration. This breaking of self-consistency is 
reminiscent of the breaking of unicity accompanying the appearence of mechanical 



instabilities for instance in elastic-plastic transition [24 
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Ref. has explored a variety of avalanches and epidemic models which have the 
same type of stress enhancement transfer at the crack tips. In noiseless systems, they 
find, in agreement with us, periodic behavior and in general large events of the size 
of the system. 



2 Description of the "crack" model 

The model is a direct extension to the crack case of the dislocation model developed in 
[llO| , m , [12| . We consider an elastic plate embedded in the (Ox, Oy) plane and composed 
of plaquettes of unit sizes paving the plane. The boundaries between the plaquettes 
constitute the elementary fault segments. They are tilted at 45 degrees with respect 
to the Ox axis, ensuring a symmetric role for all plaquette borders. A constant 
velocity boundary condition in the z direction is applied along the upper edge while 
the bottom edge is kept fixed (both the upper and bottom edges are parallel to Ox). 
Due to this externally imposed deformation and the stress transfer due to elasticity, 
each plaquette will deform. Discretizing the mechanical problem, we attribute a 
single vertical displacement w{x^ y) along the direction Oz perpendicular to the plate, 
at the center or node (x, y) of a plaquette. Each plaquette border is characterized 
by an elastic constant g which may vary from element to element (quenched disorder 
on the elastic coupling coefficients). Only two components of stress are non-zero in 
this antiplane model, namely the stress cryz{x,y) along z applied on the border/fault 
between the plaquette centered on (x, y) and the plaquette centered on (x, ?/ — 1) given 
by cTyzix, y) = g[w{x, y) — w{x,y — 1)] and the stress (Jxz{x, y) along z applied on the 
border between the plaquette centered on (x, y) and the plaquette centered on (x — 
1, y) given by cTsz(x, y) = g[w{x, y) —w{x— 1, y)]. Note that these expressions are just 
the discretized version of Hooke's law for elasticity, expressed for principal axis along 
Ox and Oy. For the present 45 degrees tilted lattice, the formulas are deduced from 
those above by the standard rule of transformation under a rigid rotation. The elastic 
displacement w{x,y) in the direction z normal to the lattice plane is solution of the 
discretized version of the equilibrium elasticity equation div(^g{x, y)gva.dw{x, y)^ = 0. 
Rupture occurs on a boundary between two plaquettes when the stress applied on it 
reaches a threshold dc which may depend on the position (quenched disorder on the 
rupture thresholds) [^. When an element breaks, the elastic strain in the element 
is relaxed but the broken element suffers no change in its material properties and 
it can support stress again in the future. The stress field is assumed to obey the 
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equation of mechanical equilibrium immediately after the rupture of an element. 
The redistribution of elastic stresses can bring other elements to rupture in a domino 
effect, creating model earthquakes. What we denominate as "fast time" is thus the 
succession of clement ruptures within an event, during which the macroscopic load at 
the plate border does not increase ("slow time" is quenched during "fast time"). This 
separation into these two time scales is intended to represent the difference between 
the fast dynamical rupture which lasts minutes at most compared to the tectonic 
loading which does not change over this time scale. 

In our previous dislocation model, the nature of the rupture on an element was 
simply characterized by the amplitude of the slip, chosen to be proportional to elastic 
deformation with a constant of proportionality f3. This amounts to model a ruptured 
element as equivalent to a dipole (antiplane is scalar) whose strength is fixed until 
the element breaks again. The total slip occurring on a fault corresponds to the 
cumulative dipole amplitude on that fault. A large earthquake in the dislocation 
model can be viewed as a nonlinear rupture pulse propagating in and being multiply 
scattered by the heterogeneous medium, with the possibility for an element to rupture 
several times in fast time. In the present "crack" version, the strength of the dipole 
is not fixed in "fast time", but must be reajusted at each rupture event in fast time 
during a given earthquake such that the dynamical stress (Tdym defined as Udyn — 
(1 — 0)(Tc where (5ac is the average stress drop, remains constant and equal to a 
preassigned value on all rupture elements in this earthquake. If n elements have 
ruptured and a new element is brought to rupture in fast time due to the stress 
redistribution induced by these n previous ruptures, the dipole strengths of the n + 1 
elements arc determined from a set of n + 1 equations as follows. Each dipole exerts a 
contribution to the stress on all other elements. The stress on any element is therefore 
the sum of the background stress prior to the earthquake plus the contribution of all 
the dipoles created in the event. These dipoles are then self-consistently determined 
such that the stress on all the ruptured elements in fast time is fixed and equal to 
the preassigned value. Physically, this models a situation where the faults remains 
"open" during the whole duration of the earthquake, at the opposite of the dislocation 
model which corresponds to an instantaneous healing (closing) of the fault after each 
rupture. 

In the crack model, there is another subtlety that was not present in the disloca- 
tion model. Suppose two or more elements are brought above their threshold in fast 
time due to the stress redistribution. The correct physics would be to solve the elas- 
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todynamics equations which direct the evolution of these unstable elements. This is 
however too difficult to implement practically for an heterogeneous system with many 
interacting faults. Our quasi-static approach circumvents this difficulty at the price 
that one has to choose, rather arbitrarily, a rule for the evolution of the unstability. 
A priori, two rules can be introduced: 1) one breaks them all simultaneously or 2) 
one ruptures only the element with the largest ratio (larger than one) of its stress 
to its threshold. We have checked that these two rules do not make any significant 
difference in the dislocation model. In the crack model, only the second rule has been 
explored in details in irreversible models of rupture, whereas the ffist rule may lead 
to un-ending ruptures. Our simulations have thus been carried out with this second 
rule. The elastic equations have been solved using a conjugate gradient technique 
with stopping criterion 10~^°. 

Most of our study will be carried out in the presence of quenched disorder in 
the stress thresholds cTc, which are drawn once for all from a probability distribution 
-P(t(c"c) chosen uniform in the interval [1 — ^,1 + ^] with the value of Act between 
0.1 and 1.9. 

The basic difference between the dislocation and crack model is that stress en- 
hancement at the fault tip is much stronger in the latter, with a stress growing as 
the square root of the fault length. As a consequence, the nature of fault and slip 
organization depends on the system size L. An earthquake of size L will generate a 
stress enhancement at its tip of magnitude equal to /SacVL- Two cases appear: 1) 
if iSacVL < A(T, the amplitude of stress enhancement generated by the dynamical 
evolution of the model is smaller than the quenched heterogeneity. The latter thus 
dominates and we expect an organization similar to that observed in the dislocation 
model where stress enhancement is small. On the contrary, for jSacVL > Aa, suffi- 
ciently large earthquakes will always create stress enhancements larger than the pre- 
existing barriers. Beyond a characteristic nucleation size L* given by (3(7cVl* ~ Act, 
earthquakes will not be stopped and will always break through the system (so-called 
"run-away" events). In addition, these large earthquakes will tend to smooth out the 
stress heterogeneity along the fault due to the condition of equal dynamical stress 
drop on the ruptured elements. These ingredients favor an approximate periodic state 
characterized by a repetition of large similar earthquakes [|TB|, 11^ EO, |23[. Take 



for instance (3 = 0.1. Then, the run-away will be absent for systems sizes smaller 
than of the order of 100. For such small stress drops, we have verified in particular 
that the fault patterns selected in the dislocation and the crack models have similar 
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statistical properties in this regime. For larger stress drops or larger system sizes, 
there are two populations of events : 1) the small ones similar to that occurring in 
the dislocation model albeit with a different distribution; 2) large earthquakes of size 
comparable to the system size. 

3 Threshold disorder 

Here, we wish to explore the crack regime. Our simulations have thus been carried 
out on large systems 130 by 130. The influence of the dynamical stress a^yn has been 
explored in the range < adyn < I — since otherwise unending ruptures occur 
(the second inequality expresses that the dynamical stress must obviously be smaller 
than all static stress thresholds). 

• Aa = 0.1 (leading to adyn < 0.95). 

For adyn = (large stress drop), after a short transient, the dynamics is that of a 
perfect "characteristic earthquake" |^: because of our tilted lattice structure, a well- 
deflned regular fault, made of two linear strands oriented at 45 degrees with respect to 
the Ox axis and forming a V, is activated regularly in a perfect periodic fashion by a 
single great earthquake in which all elements on the fault break once in fast time with 
exactly the same slip. This regime corresponds to a perfect synchronization of all the 
threshold elements constituting the fault. The Gutenberg- Richter distribution is a 
Dirac function. The same is found for intermediate stress drop {adyn = 0.4). For small 
stress drop {adyn = 0.9), we flnd again a perfect synchronization corresponding to the 
repetition of a single large identical event. However, the fault on which this event 
occurs is now rough, characterized by linear strands separated by rough portions. Its 
speciflc structure is however dependent upon the speciflc realization of the disorder 
on the thresholds. In summary, a small disorder favors a very regular organization. 

• Aa = 1 (leading to adyn < 0.5). For adyn = 0, the behavior is very similar 
to the previous case Aa = 0.1, except that the fault is now rough all along its 
length. For larger dynamical stress adyn < 0.2 and 0.4, the dynamics is still periodic. 
However, a period contains a much more complex history than just the succession, 
as documented until now, of a quiescent phase followed by a single great earthquake. 
In contrast, after a great earthquake, there is long quiescent phase, followed by the 
appearence of a diffuse foreshock activity spread over the plate. This diffuse activity 
is made of many small earthquakes. It accelerates up to the time where the great 
earthquake occurs on a fault. This fault is again well-defined and is finally selected 
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after a long transient. The distribution of earthquake energies contains two parts: a 
nice powerlaw distribution P{E)dE ~ for small earthquakes 0.2 < £^ < 20 

(in the units where the elastic coefficients are all equal to 1) with B = 1 ± 0.05 and 
a Dirac peak at the energy of the great event (around E ~ 4000). Notice the huge 
separation of energy scales. This can be rationalized using the nucleation argument 
outhned above. From the expression PacVL* ~ Act, we get a nucleation scale L* of 
the order of 4 elements. The energy being proportional to the square of the length in 
crack elasticity, this yields a characteristic maximum energy scale of 16 not far from 
the maximum energy observed in the power law distribution. In contrast, the great 
earthquake has a size of the order of one hundred and its release energy is thus of 
order 10^. 

• Act = 1.9 (leading to a^yn < 0.05). The progressive organization of the earth- 
quake activity on a localized fault network is shown in fig.l. This case is similar 
to the previous one with the larger a^yn- However, quantitatively, the rupture his- 
tory during a period separating the recurrence of two great earthquakes is even more 
characterized by the appearence of a multitude of small earthquakes (see fig. 2c). The 
activity is always localized on a well-defined fault structure (see fig. 2a) which becomes 
fixed at large times. However, the fault is no more a linear object, but contains loops 
defining internal " micro-plates" . For the largest allowed dynamical stress drops (say 
(^dyn — 0.04), we observe additional properties. After a great earthquake, there is, 
as usual, a quiescent time, followed by a progressing activation of the main fault by 
small earthquakes (fig.2c). The main fault is defined as the locus of the great earth- 
quakes (see fig.2a). The activity of small earthquakes then shifts progressively away 
from the main fault to become delocalized in the bulk of the plate. This activity 
accelerates until the great earthquake occurs on the main fault. Again, we observe a 
huge separation of energy scales between the small events distributed according to a 
power law distribution for energies between 10~^ and a few 10~^ and the great char- 
acteristic earthquake of energy of the order of 1500. This is again correctly explained 
by the nucleation argument. The exponent of the power law distribution for small 
earthquakes is again B — 1± 0.05 (see fig.2d) 

Before the periodic regime organizes itself (see fig. 2b) with its complex spatio- 
temporal structure of small earthquakes preceding the run-away, we witness a long 
transient aperiodic regime. The time duration of this transient is all the longer, the 
larger is the dynamical stress (Tdyn- For adyn = 0.04 for instance, the transient is so 
long that one can measure with good accuracy the distribution of earthquake energies 
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in time windows sufficiently large so that the statistics is good, but sufficiently small 
so that the evolution of the organization in this transient regime is negligible. (Note 
however that the transient is nevertheless quite small compared to that observed 
in the dislocation model: the crack model introduces a much stronger coupling via 
the crack stress redistribution rule, favoring a faster and stronger organization). We 
again obtain a nice power law for small earthquakes in the transient regime, with the 
same exponent B. However, the sizes of these "small" earthquakes in the transient 
regime are typically one order of magnitude larger than in the asymptotic periodic 
regime. 

In summary, these numerical explorations show that the asymptotic dynamics is 
always periodic, with however a more and more complex sequence of ruptures within 
a period, the larger are the disorder Act and the dynamical stress adyn (i-e. the 
smaller the stress drop). This can be qualitatively understood as an intermediate 
regime between the fully periodic regime, which is characteristic of low disorder and 
large stress drop, and the self-organized critical behavior observed for the dislocation 
model, which is controlled essentially by large disorder and small stress drop [lin. 



ni IT2[. We have already underlined the analogy between this problem and that 



of a set of interacting relaxation oscillators with threshold. In this analogy, the 
stress drop parameter measures the coupling strength between elements, whereas the 
amplitude of the disorder Act quantifies the disorder in the natural frequencies of 
the individual elements. The analogy predicts that synchronization, hence regular 
periodic behavior, will be the stronger the stronger the coupling and the weaker the 



disorder 18, m m 21, 22, E 



Let us present a mean field toy version of the model, in the spirit of [Q, which 
allows one to understand the mechanism underlying the synchronization process and 
the appearence of periodic states. For the purpose of clarity, consider an homogeneous 
system and the situation where the earthquake cycle is constituted of two earthquakes, 
involving respectively A^i and N2 < Ni fault elements. The mean field character of the 
argument is to assume that, when an element reaches its threshold Uc, it redistributes 
its stress to all the other N = Ni + N2 active elements, each thus getting a stress 
increment equal to where we allow for an arbitrary positive coupling strength a. 
For simplicity of the argument, suppose also that the stress of the ruptured element 
is put to zero. Suppose that we start reasoning at the time where the stress on the 
elements of fault 1 is di > (J2, where (T2 is the stress on the elements of fault 2. 
Earthquake 1 will occur ffist, when all A^^i elements are at their ac in fast time. At 



10 



that time, the stress on the elements of fault 2 is cr2 + o^c — en, due to the uniform 
tectonic loading. Just after event 1, the stress on the elements of 1 is zero by our 
rules while the stress on the fault 2 is o"2 + — o"i + . Fault 2 will then reach 



N 

itself (Tc at which time the stress on fault 1 is ai — a2 — After earthquake 2, 

the stress on fault 2 is zero by definition while the stress transfer loads fault 1 to the 
level ai — (12 — iEi:zI^i}£!£^ _ This last result shows that the stress difference, which was 
initially ai — has decreased during an earthquake cycle. The largest earthquake 
(1 in this example) is an absorbing state. This argument summarized here for an 
homogeneous system works also for heterogeneous systems, if the disorder is not too 
large . 



4 Elasticity and threshold disorder 

We explore in this section the possibility to destroy the periodic behavior and kill 
the great earthquake, by introducing disorder also in the elastic coefficients of the 
elements and by varying the nature of the disorder, for instance by allowing for the 
existence of very strong and rigid elements (described by a powerlaw distribution of 
thresholds and/or elastic coefficients). This comes about because, the stronger the 
disorder, the stronger will be the barriers to stop the run-away. 

With respect to the fault structure and the time history, the addition of disorder 
on the elastic coefficient, when not too large, is tantamount to increasing the threshold 
disorder with no elastic disorder: we observe a periodic behavior, after a transient 
which is significantly larger than previously, all parameters being the same otherwise. 
When measurable, the distribution of small earthquakes is a power law with an 
exponent B always equal to 1±0.05. Fig. 3 presents such a simulation with Act = 1.9, 
a disorder on the elastic coefficients defined as for the threshold with a flat distribution 
of elastic coefficient and a width equal to l\g = 1 and adyn = 0.04. 

For large threshold disorder Act = 1.9, we find the novel feature that, in some 
system realizations, the time history does not seem to become periodic, however long 
we wait. A global stationary regime seems to emerge nevertheless, with an elastic 
energy stored in the plate which fluctuates around a well-deflned average. The fault 
structure becomes very rough with overhangs (this is allowed in the scalar anti-plane 
elastic model used here but would be unphysical in in-plane stress or strain models 
since compressive and extensive stress would accumulate without limit in the region 
of the overhang, according to a "hook" effect 0]. Furthermore, the run-away does 
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not exist anymore and is replaced by a continuous distribution of earthquake of all 
sizes. 

We have also explored different disorders, for instance a powerlaw distribution 
of rupture thresholds Psigma{(^c) ~ cr^(i+'^'^) for dc > 1, with /Xg. = 0.5 and 3. The 
first case corresponds to a very broad distribution (the average is mathematically 
divergent) with a significant fraction of the bonds having a very large threshold. De- 
pending up the specific realization, we find a similar phenomenology as before. For 
instance, for adyn = 0, a single great run-away punctuated the periodic dynamics. 
However, the fault becomes very complex, with many branches and loops. Correla- 
tively, the number of elements rupturing in the great earthquake is about three times 
the system width. 

For other disorder realizations with the powerlaw distribution of rupture thresh- 
olds, in presence or absence of elastic disorder, we have found a completely new 
effect, that we identify as a "frustration" which destroys the self-consistency of the 
quasi-static crack model. This also occurs for the previous bounded distribution in 
the presence of quenched disorder in the elastic coefficients. The phenomenon is best 
illustration by examination of figure 4. It shows a small part of the stress field during 
an earthquake in fast time in a network with a power law distribution of thresholds 
with /io- = 0.5 and adyn = 0.5. The fault segments which have ruptured are indicated 
by an arrow with a thick line. These ruptured elements all carry a stress whose 
absolute value is (Tdyn = 0.5. The arrows indicate the sign of the stress: downward 
(resp. upward) arrows correspond to positive (resp. negative) stresses (recall that 



all the stresses are along the Oz axis) p9[. The stress carried by each fault element 
is indicated on the arrow. The rupture threshold for each element is written on the 
other side of the bond. The mechanical equilibrium translates into the condition that 
the sum of stresses carried by the fault segments attached to the same node be zero, 
easily verified in this example. The self-consistency of the model is obeyed when the 
mechanical equilibrium is such that all stresses are smaller than the corresponding 
threshold. When this is not the rupture occurs relaxing the stress on the 

rupture element to adyn- However, there are "frustrated" situations in which this is 
not possible. To see this, examine the element carrying a stress equal to 1.5 whose 
threshold is cTc = 1.495. By the definition of the model, since its stress is larger than 
its threshold, it must rupture and its stress must thus be lowered to adyn = 0.5. When 
this occurs, the mechanical equilibrium is broken, because the three other fault ele- 
ments connected to the same node have all their stresses imposed to be equal to the 
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dynamical stress. The model looses its self-consistency due to an over-determination. 
In general, it is straightforward to realize that this "frustrated" state occurs when 
three fault elements connected to the same node have ruptured. As a consequence, 
their stresses is imposed to be equal to adyn and from the rule of quasi-static equilib- 
rium, the stress on the remaining element connected to the same node is completely 
determined. If it happens that its threshold be less than this stress, the frustration 
appears and the model collapses. 

This frustration is reminiscent of the concept of frustration introduced in the 



physics of general disordered systems |31[]. In general, frustration arises in a system 
whose interactions compete or conflict in such a way that not all constraints on the 
system can be simultaneously satisfied. This is what happens in our model with 
the conflict between the stress conservation and the dynamical stress which may 
produce frustration. The usual outcome of frustation in the presence of disorder is 
the existence of, not a single well-defined stable equilibrium state, but rather to a 



very large number of disordered equilibrium states of equivalent energies |^2[. The 
number of these minimum states usually increases exponentially with the number 
of degrees of freedom (spins). The energy landscape is extremely complicated with 
a hierarchy of barriers of increasing sizes separating the minimum states. In other 
words, in order to go from one minimum state to another, an energy barrier must be 
passed. This hierarchical structure produces novel behaviors (long time relaxations, 
breakdown of ergodicity, etc.) and important fluctuations. In our model, mechanical 
equilibrium is not more possible in the presence of frustration, as we have shown. It 
is tempting to extrapole from this analogy and suggest that the state of stress chosen 
by the plate, in an extended model allowing for a genuine dynamics, could have a 
multi- valley structure similar to that obtained in spinglasses for instance. This could 
be an important underlying ingredient at the origin of the observed spatio-temporal 
complexity. In this spirit, we have recently proposed that the mechanics of coupled 
blocks and especially the coupling between rotations exhibit the frustration property 
|]33| . We believe that this is an ubiquitous and key property at the basis of the 
complexity observed in fault and earthquake organization. 

The present crack model must be enriched to get rid of the overdetermination. 
The conceptually simplest solution is to introduce a genuine dynamics allowing a 
transient breakdown of static equilibrium. For instance, elastodynamics allow for 
an unbalance of the stress, corresponding to the generation and radiation of elastic 
waves. Physically, this describes the fact that the stored elastic energy is converted 
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into two types of dissipation: 1) frictional heating, that we take into account by our 
rupture criterion and 2) elastic wave radiation that we neglect. A model in which the 
dissipation is not completely converted into friction, i.e. the stress is not immediately 
put to its asymptotic dynamical value, would cure the observed frustration and over- 
determination. Another alternative, within quasi-static mechanics, is to introduce 
other types of deformations, such as thrust or normal faulting. Algorithmically, this 
will lead to allow a fraction of the stress to be removed by going out of plane, thereby 
removing the constraint of zero antiplane stress at all nodes. However, we still expect 
frustation to occur for certain realizations due to the competition between static 
mechanical equilibrium and constant dynamical stress. 

It is interesting to note that there is a particular value of the dynamical stress drop 
for which the frustration never occurs, namely for adyn = 0. In this case, the stress 
on the fourth element attached to a node, for which the three other fault segments 
have ruptured, is zero by the law of local static mechanical equilibrium. This is 
consistent both with mechanical equilibrium by definition and also with the dynamical 
stress drop condition. Notice that this is the only situation for which the condition 
of dynamical stress drop is always compatible with static mechanical equilibrium. 
However, the crack model presents a curious and probably rather artificial behavior 
in this case. Remember that a rupture cycle is characterized usually by a progressive 
acceleration of the rate of small earthquakes prior to the occurrence of the run-away. 
Since the run-away spans the whole width of the lattice, and since the dynamical 
stress is imposed on its ruptured elements, this amounts to impose the total stress 
within the plate equal to zero. Previous ruptures on small faults within a cycle have 
produced localized slip and stress sources. They cannot remain in the presence of 
this global vanishing of the stress within the plate. As a consequence, we observe 
in fast time a cloud of small earthquakes accompanying the run-away, which are the 
"ghosts" of all the previous small earthquakes. These "ghosts" present exactly a slip 
which is the opposite of the shp that they have develop in the foreshock phase. In 
other words, the aftershocks occurring in fast time are the exact symmetric of all the 
foreshocks. The difference is that the foreshocks are spread in time over the period of 
the cycle while the aftershocks "ghosts" are occurring in fast time, just after the run- 
away. While the details of this behavior is clearly model specific, this phenomenon 
is not without recalling field observations that foreshocks occur usually years or even 
decades before a great event, while the huge majority of aftershocks are clustered over 
a few months after the main event. The present model does not contain however the 
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necessary ingredients to describe the time delays associated with the couphng with 
other modes of creep or ductile deformations, the ductile crust and the fluid in the 
crust. 



Ref. |3J] have studied the same crack model and it is instructive to discuss how 
their results differs from ours in many aspects. They have only studied the case 
<^dyn = 0. They have thus not found the frustration effect discussed above. In 
addition, they consider an annealed disorder, i.e. all the threshold of the fractures 
elements are re-set to new random numbers after each event. As a consequence, they 
cannot obtain earthquake localization on well-defined faults but only observe diffuse 
earthquakes. Also, this reshuffling of the disorder prevents the synchronization to a 
periodic cycle and the appearence of a run-away. Nevertheless, they observe that the 
distribution of small earthquakes is a powerlaw with B = 0.8 ± 0.1, not far from our 
estimate and that the distribution presents a peak at large earthquakes whose energy 
scales with the square of the system size. While they interpret this as a finite size 
effect, we rather conclude that these large events are the shadows of the great run- 
aways in the presence of annealed noise, which sizes are given by that of the system. 
In a sense, the annealed disorder makes their system function permanently in our 
transient regime. They have only studied disorder on the thresholds by the method 
of Green functions, which is not useful practically in the presence of elastic disorder. 
The gradient conjugate method that we have used is slower but more general to 
tackle this second case. Finally, they have used infinite system Green functions, and 
have not addressed the question of the effect of boundaries in finite systems. While 
in the statistical physics of critical phenomena, one would like to get results which 
are independent of boundaries, in the present mechanical problem as well as in the 
general mechanical case, the existence of well-defined boundary conditions on the 
stress or strain fields at the border of the system is know to control drastically the 
localization of the mechanical deformation, as we have been able to observe. 



5 Concluding remarks 

Except for special realizations with the strongest disorder on rupture thresholds and 
elastic coefficients, we have found that the quasi-static crack model of earthquake 
recurrence leads to periodic cycles, characterized by small foreshocks distributed ac- 
cording to a universal Gutenberg-Richter law with exponent i? = 1 up to a maximum 
nucleation size and a large run-away ending the cycle. This periodic behavior results 
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from the strong synchronization brought by imposing a constant dynamical stress, 
corresponding to an attractive or absorptive state. The other main result is the 
discovery of a fundamental frustration resulting from an overdetermination of the 
stress field in the presence of large disorder and imposed dynamical stress drop. The 
general solution to this breakdown of self-consistency is to re-introduce a genuine 
dynamics allowing the local breakdown of static mechanical equilibrium, associated 
to the radiation of elastic waves. Our study pinpoints the fundamental role played 
by elastodynamics in repetitive crack ruptures. We thus believe that, in crack mod- 
els (and not in dislocation models), there is no other way than incorporate the full 
elastodynamic equations to get a self-consistent solution in all situations. 
We acknowledge stimulating discussions with L. Knopoff and S. Nielsen. 
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Figure captions 

Fig.l : The parameters are Aa = 1.9, a^yn = 0.04, and the constant velocity 
imposed at the boundary isV = 10^^. The figure shows two maps of the accumulated 
slip in fault segment at two different times in a square system of size L = 130 in the 
transient regime. We represent those elements which have slipped at least once, the 
light grey to black scale corresponding to increasing cumulative slips. The two times 
of observations correspond respectively to the first 1500 events (top) and 2000 events 
(bottom) . One observes a rather diffuse " damage" at early times and localization of 
the deformation at long times. 

Fig. 2 : The parameters are Act = 1.9, a^yn — 0.04, and the constant velocity 
imposed at the boundary isV— 10~^. 

a) Long-time accumulated slip after localization for the same parameters as for 
fig.l but with another disorder reahzation. All the earthquakes occur on this ID 
fault. 

b) Time dependence of the total elastic energy stored in the plate. The periodic 
behavior established after the transient up to time 50000 is clearly visible. 

c) Coding of the active fault elements as a function of time. The position of a 
given fault element is coded by a single index (denoted "position of broken bonds" 
in the figure) increasing from 1 to = 16900. The index spans x = 1 to L at fixed 
y for each y = 1 to L. Note the periodic cycle with three phases: i) quiescence after 
the great event, ii) reactivation of the seismic activity mainly on the main fault and 
iii) increase of the foreshocks frequency in the system. 

d) Gutenberg-Richter distribution of the number of events having a given energy. 
The energy of an event is defined as the difference between the total elastic energy 
stored in the system before and after the event. The crosses, plus and square symbols 
corresponds to different times in the dynamics at long times, showing the stability of 
the distribution. The diamonds correspond to the distribution of small events in a 
time window in the third regime close to the run-away in the periodic regime. In all 
cases, we observe powerlaw distribution P{E)dE ~ with a constant exponent 
B — 1 :k 0.05. The straight fine has slope —2 for comparison. It is interesting to 
note that the foreshocks represented by the diamonds have the same S- value but are 
larger on average, signaling the nucleation of the run-away. The great earthquake is 
not represented on the figure, being out of scale. 
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Fig. 3 : a) Fault localization obtained at long times for Au = 1.9, disorder on the 
elastic coefficients = 1, aayn = 0.04, and the constant velocity imposed at the 
boundary is V = 10~^. Note that the lower fault in grey becomes inactive after the 
transient and its contribution to the cumulative slip vanishes at long times. 

b) Gutenberg-Richter distribution of the number of events having a given energy 
for different time windows at increasing times when going from right to left. In all 
cases, we observe powerlaw distribution P{E)dE ~ with a constant exponent 

B — 1± 0.05. The straight line has slope —2 for comparison. For early times, in 
the power law presents the same exponent, the events have a larger size which finally 
settle to a stationary distribution. 

Fig.4 : Map of a small part of the stress field during an earthquake in fast time in 
a network with a power law distribution of thresholds with //g. = 0.5 and adyn = 0.5. 
The fault segments which have ruptured are indicated by an arrow with a thick line. 

These ruptured elements all carry a stress whose absolute value is adyn = 0.5. The 
arrows indicate the sign of the stress: downward (resp. upward) arrows correspond to 
positive (resp. negative) stresses (recall that all the stresses are along the Oz axis). 
The stress carried by each fault element is indicated on the arrow. The rupture 
threshold for each element is written on the other side of the bond. A "frustrated" 
element is seen, with a stress equal to 1.5 whose threshold is (Tc = 1.495. 
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